
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> 
> # change directory to the replication location
> setwd("D:/Dropbox/jesarey_documents/PS Article on Reviewing")
> 
> ###########
> # combined plot
> # unanimity, majority, majority w/ editor, and unilateral editor rules
> # 10% reviewer/editor acceptance rate
> 
> # 3 reviewers, unanimity required
> 
> rm(list=ls())
> set.seed(123456)
> 
> library(copula)
> 
> rho <- seq(from=0, to=0.98, by=0.02)
> pr.accept <- 0.1
> 
> accept<-c()
> for(k in 1:length(rho)){
+   reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k]), dim=3, dispstr="un")) >= (1-pr.accept)
+   decisions <- apply(X=reviews, MARGIN=1, FUN=min)
+   
+   # acceptance rate
+   accept[k] <- sum(decisions)/length(decisions)
+ }
> 
> png(file="acceptance-rates.png", width=5, height=6, units="in", res=144)
> plot(accept ~ rho, ylim=c(0, 0.13), col=gray(0.5), ylab = "proportion of accepted manuscripts", xlab = "correlation of reviewer opinions", las=2)
> y.fit <- predict(loess(accept~rho))
> lines(y.fit~rho, lty=3, lwd=2)
> 
> 
> 
> # 3 reviewers, majority voting
> 
> rm(list=ls())
> set.seed(123456)
> 
> rho <- seq(from=0, to=0.98, by=0.02)
> pr.accept <- 0.1
> 
> accept<-c()
> for(k in 1:length(rho)){
+   reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k]), dim=3, dispstr="un")) >= (1-pr.accept)
+   decisions <- ifelse( apply(X=reviews, MARGIN=1, FUN=sum) >= 2, 1, 0 )
+   
+   # acceptance rate
+   accept[k] <- sum(decisions)/length(decisions)
+ }
> 
> points(accept ~ rho, pch=5, col=gray(0.5))
> y.fit <- predict(loess(accept~rho))
> lines(y.fit~rho, lwd=2, col=gray(0.5))
> 
> 
> # editor as fourth reviewer, 3 out of 4 positive reviews required
> 
> rm(list=ls())
> set.seed(123456)
> 
> rho <- seq(from=0, to=0.98, by=0.02)
> pr.accept <- 0.1
> 
> accept<-c()
> for(k in 1:length(rho)){
+   reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k], rho[k], rho[k], rho[k]), dim=4, dispstr="un")) >= (1-pr.accept)
+   decisions.t <- apply(X=reviews, MARGIN=1, FUN=sum)
+   decisions <- decisions.t>=3
+   
+   # acceptance rate
+   accept[k] <- sum(decisions)/length(decisions)
+ }
> 
> points(accept ~ rho, pch=4, col=gray(0.5))
> y.fit <- predict(loess(accept~rho))
> lines(y.fit~rho, lty=2, lwd=2)
> 
> 
> # strong editor system with truthful reports
> 
> rm(list=ls())
> set.seed(123456)
> 
> rho <- seq(from=0, to=0.98, by=0.02)
> pr.accept <- 0.1
> 
> accept<-c()
> for(k in 1:length(rho)){
+   reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k], rho[k], rho[k], rho[k]), dim=4, dispstr="un"))
+   decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
+   decisions <- decisions.t >= (1-pr.accept)
+   
+   # acceptance rate
+   accept[k] <- sum(decisions)/length(decisions)
+ }
> 
> points(accept ~ rho, pch=2, col=gray(0.5))
> y.fit <- predict(loess(accept~rho))
> lines(y.fit~rho, lwd=1)
> 
> abline(h=0, col="lightgray")
> 
> source(file="LEGEND.r")
> 
> # source: http://stackoverflow.com/questions/19053440/r-legend-with-points-and-lines-being-different-colors-for-the-same-legend-item
> LEGEND("topleft", bty="n", lty=c(1,2,1,3), pch=c(2,4,5,1), pt.lwd=1, cex=0.8,  lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col="black", line.col=c("black", "black", gray(0.5), "black"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> # 
> # Very large reviewer pool
> ##########################
> 
> rm(list=ls())
> library(copula)
> set.seed(123456)
> rho <- 0.5
> 
> accept<-c()
> 
> reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
> reviews <- reviews.all[,1:4]
> 
> # strong editor system
> decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.2)
> sum(decisions)/length(decisions)
[1] 0.1058
> 
> # unanimity voting
> decisions.t.2 <- reviews[,1:3] >= (1-0.3)
> decisions.2 <- apply(X=decisions.t.2, MARGIN=1, FUN=min)
> sum(decisions.2)/length(decisions.2)
[1] 0.10012
> 
> # majority voting with active editor
> decisions.t.3 <- reviews >= (1-0.2)
> decisions.3 <- apply(X=decisions.t.3, MARGIN=1, FUN=sum) >=3
> sum(decisions.3)/length(decisions.3)
[1] 0.10052
> 
> # majority voting
> decisions.t.4 <- reviews[,1:3] >= (1-0.15)
> decisions.4 <- apply(X=decisions.t.4, MARGIN=1, FUN=sum) >=2
> sum(decisions.4)/length(decisions.4)
[1] 0.11112
> 
> reviews.avg <- apply(reviews.all, MARGIN=1, FUN=mean)
> 
> png(file="opinion-density.png", width=3.6, height=4.3, units="in", res=144)
> 
> plot(density(reviews.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,6), xlab="mean reader quality evaluation (percentile)", main="Panel 2A: Without Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
> lines(density(reviews.avg[which(decisions.3 == T)], from=0, to=1), lwd=2, lty=2)
> lines(density(reviews.avg[which(decisions.4 == T)], from=0, to=1), lwd=3, col=gray(0.5))
> lines(density(reviews.avg[which(decisions.2 == T)], from=0, to=1), lwd=2, lty=3)
> 
> legend("topleft", bty="n", lty=c(1,2,1,3), cex=0.6, lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col=c("black", "black", gray(0.5), "black"))
> dev.off()
null device 
          1 
> 
> # proportion of majority voting results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.4 == T)])(0.65)
[1] 0.1171706
> # proportion of majority voting with active editor results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.3 == T)])(0.65)
[1] 0.06884202
> # proportion of unanimity voting results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.2 == T)])(0.65)
[1] 0.1132641
> # proportion of strong editor results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions == T)])(0.65)
[1] 0.0557656
> 
> mean(reviews.avg[which(decisions.4==T)])
[1] 0.7767603
> mean(reviews.avg[which(decisions.3==T)])
[1] 0.7943731
> mean(reviews.avg[which(decisions.2==T)])
[1] 0.779667
> mean(reviews.avg[which(decisions==T)])
[1] 0.7976532
> 
> # proportion of majority voting results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.4 == T)])(0.90)
[1] 0.8929086
> # proportion of majority voting with active editor results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.3 == T)])(0.90)
[1] 0.879228
> # proportion of unanimity voting results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.2 == T)])(0.90)
[1] 0.8893328
> # proportion of strong editor results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions == T)])(0.90)
[1] 0.8833648
> 
> png(file="acceptance-prob.png", width=3.6, height=4.3, units="in", res=144)
> 
> # make a plot showing the role of "chance" in the review process
> avg.quality <- apply(X=reviews.all, MARGIN=1, FUN=mean)
> prob.fit <- loess(decisions.4 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> plot(prob.accept ~ avg.qual.fit, lwd=2, col=gray(0.5), type="l", ylim=c(0, 1), xlab = "mean reader quality evaluation (percentile)", ylab = "probability of acceptance", main="Panel 3A: Without Desk Rejection", yaxs="i", cex.main=0.9, las=2, cex.lab=0.9)
> 
> 
> prob.fit <- loess(decisions.2 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit, lty=3, lwd=2)
> 
> 
> prob.fit <- loess(decisions.3 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit, lty=2, lwd=2)
> 
> 
> prob.fit <- loess(decisions ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit)
> 
> legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice        ", "majority rule, reviewers + editor        ", "majority rule, reviewers only        ", "unanimity rule, reviewers only         "), col=c("black", "black", gray(0.5), "black"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> # impose desk rejections for the papers the editor thinks are the worst
> not.desk <- which(reviews[,4] >= 0.5)
> 
> # strong editor system
> decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.13)
> sum(decisions)/length(decisions)
[1] 0.09645056
> 
> # unanimity voting
> decisions.t.2 <- reviews[not.desk,1:3] >= (1-0.21)
> decisions.2 <- apply(X=decisions.t.2, MARGIN=1, FUN=min)
> sum(decisions.2)/length(decisions.2)
[1] 0.09956831
> 
> # majority voting with active editor
> decisions.t.3 <- reviews[not.desk,] >= (1-0.15)
> decisions.3 <- apply(X=decisions.t.3, MARGIN=1, FUN=sum) >=3
> sum(decisions.3)/length(decisions.3)
[1] 0.1253098
> 
> # majority voting
> decisions.t.4 <- reviews[not.desk,1:3] >= (1-0.09)
> decisions.4 <- apply(X=decisions.t.4, MARGIN=1, FUN=sum) >=2
> sum(decisions.4)/length(decisions.4)
[1] 0.1051243
> 
> reviews.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)
> 
> png(file="opinion-density-desk.png", width=3.6, height=4.3, units="in", res=144)
> 
> plot(density(reviews.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,6), xlab="mean reader quality evaluation (percentile)", main="Panel 2B: With Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
> lines(density(reviews.avg[which(decisions.3 == T)], from=0, to=1), lwd=2, lty=2)
> lines(density(reviews.avg[which(decisions.4 == T)], from=0, to=1), lwd=3, col=gray(0.5))
> lines(density(reviews.avg[which(decisions.2 == T)], from=0, to=1), lwd=2, lty=3)
> 
> legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col=c("black", "black", gray(0.5), "black"))
> 
> dev.off()
null device 
          1 
> 
> # proportion of majority voting results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.4 == T)])(0.65)
[1] 0.03840304
> # proportion of majority voting with active editor results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.3 == T)])(0.65)
[1] 0.03508772
> # proportion of unanimity voting results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions.2 == T)])(0.65)
[1] 0.03412284
> # proportion of strong editor results at less than the 65th percentile
> ecdf(reviews.avg[which(decisions == T)])(0.65)
[1] 0.01409034
> 
> mean(reviews.avg[which(decisions.4==T)])
[1] 0.8244604
> mean(reviews.avg[which(decisions.3==T)])
[1] 0.8213632
> mean(reviews.avg[which(decisions.2==T)])
[1] 0.8242498
> mean(reviews.avg[which(decisions==T)])
[1] 0.8405842
> 
> # proportion of majority voting results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.4 == T)])(0.90)
[1] 0.8007605
> # proportion of majority voting with active editor results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.3 == T)])(0.90)
[1] 0.8223285
> # proportion of unanimity voting results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions.2 == T)])(0.90)
[1] 0.806102
> # proportion of strong editor results at less than the 90th percentile
> ecdf(reviews.avg[which(decisions == T)])(0.90)
[1] 0.7733112
> 
> png(file="acceptance-prob-desk.png", width=3.6, height=4.3, units="in", res=144)
> 
> # make a plot showing the role of "chance" in the review process
> avg.quality <- reviews.avg
> prob.fit <- loess(decisions.4 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> plot(prob.accept ~ avg.qual.fit, lwd=2, col=gray(0.5), type="l", ylim=c(0, 1), xlab = "mean reader quality evaluation (percentile)", ylab = "probability of acceptance", main="Panel 3B: With Desk Rejection", yaxs="i", cex.main=0.9, las=2, cex.lab=0.9)
> 
> 
> prob.fit <- loess(decisions.2 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit, lty=3, lwd=2)
> 
> 
> prob.fit <- loess(decisions.3 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit, lty=2, lwd=2)
> 
> 
> prob.fit <- loess(decisions ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
> avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
> prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
> lines(prob.accept ~ avg.qual.fit)
> 
> legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice        ", "majority rule, reviewers + editor        ", "majority rule, reviewers only        ", "unanimity rule, reviewers only         "), col=c("black", "black", gray(0.5), "black"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> ###########
> # plot average review quality against percentile
> # strong editor system with truthful reports
> # two subfield discipline
> 
> rm(list=ls())
> set.seed(123456)
> 
> library(copula)
> 
> N <- 250              # size of each subdiscipline
> 
> rho.hi <- 0.9         # within-subfield correlation
> rho.lo <- 0.1         # between-subfield correlation
> 
> rho.vec <- c()
> for(i in (N-1):0){
+   rho.vec <- c(rho.vec, rep(rho.hi, i), rep(rho.lo, N))
+ }
> rho.vec <- c(rho.vec, rep(rho.hi, sum(1:(N-1))))
> 
> reviews.all <- rCopula(50000, normalCopula(param=rho.vec, dim=(N*2), dispstr="un"))
> # calculate average correlation
> cor.mat <- cor(reviews.all)
> mean(cor.mat)
[1] 0.4962566
> 
> reviews <- reviews.all[,c(1:2, 499:500)]
> 
> # strong editor system
> decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.22)
> sum(decisions)/length(decisions)
[1] 0.1016
> 
> review.avg <- apply(X=reviews, MARGIN=1, FUN=mean)
> rank.quant <- ecdf(review.avg)(review.avg)
> o<-order(rank.quant)
> 
> png(file="opinion-density-subdis.png", width=3.6, height=4.3, units="in", res=144)
> 
> plot(density(review.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,13.5), xlab="mean reader quality evaluation (percentile)", main="Panel 4A: Without Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.8536645
> 
> 
> # repeat with one field only, low correlation
> rm(list=ls())
> set.seed(123456)
> rho <- 0.5
> 
> accept<-c()
> 
> reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
> reviews <- reviews.all[,1:4]
> 
> # strong editor system
> decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.2)
> sum(decisions)/length(decisions)
[1] 0.1058
> 
> review.avg <- apply(reviews.all, MARGIN=1, FUN=mean)
> 
> lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lty=2)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.7976532
> 
> 
> 
> 
> # repeat with one field only, high correlation
> rm(list=ls())
> set.seed(123456)
> rho <- 0.75
> 
> accept<-c()
> 
> reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
> reviews <- reviews.all[,1:4]
> 
> # strong editor system
> decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.15)
> sum(decisions)/length(decisions)
[1] 0.10498
> 
> review.avg <- apply(reviews.all, MARGIN=1, FUN=mean)
> 
> lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lwd=2)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.8847323
> 
> 
> 
> 
> legend("topleft", bty="n", cex=0.6, lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black"), legend=c("two subdisciplines, rho = 0.9 or 0.1          ", "all readers correlated at rho = 0.5          ", "all readers correlated at rho = 0.75          ") )
> 
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> ###########
> # plot average review quality against percentile
> # strong editor system with truthful reports
> # two subfield discipline
> # with desk rejection
> 
> rm(list=ls())
> set.seed(123456)
> 
> library(copula)
> 
> N <- 250              # size of each subdiscipline
> 
> rho.hi <- 0.9         # within-subfield correlation
> rho.lo <- 0.1         # between-subfield correlation
> 
> rho.vec <- c()
> for(i in (N-1):0){
+   rho.vec <- c(rho.vec, rep(rho.hi, i), rep(rho.lo, N))
+ }
> rho.vec <- c(rho.vec, rep(rho.hi, sum(1:(N-1))))
> 
> reviews.all <- rCopula(50000, normalCopula(param=rho.vec, dim=(N*2), dispstr="un"))
> reviews <- reviews.all[,c(1:2, 499:500)]
> 
> # impose desk rejection
> not.desk <- which(reviews[,4] >= 0.5)
> 
> # strong editor system
> decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.15)
> sum(decisions)/length(decisions)
[1] 0.09456369
> 
> review.avg <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
> rank.quant <- ecdf(review.avg)(review.avg)
> o<-order(rank.quant)
> 
> png(file="opinion-density-subdis-desk.png", width=3.6, height=4.3, units="in", res=144)
> 
> plot(density(review.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,13.5), xlab="mean reader quality evaluation (percentile)", main="Panel 4B: With Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.9008516
> 
> 
> # repeat with one field only, low correlation
> rm(list=ls())
> set.seed(123456)
> rho <- 0.5
> 
> accept<-c()
> 
> reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
> reviews <- reviews.all[,1:4]
> 
> # impose desk rejection
> not.desk <- which(reviews[,4] >= 0.5)
> 
> # strong editor system
> decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.13)
> sum(decisions)/length(decisions)
[1] 0.09645056
> 
> review.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)
> 
> lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lty=2)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.8405842
> 
> 
> # repeat with one field only, higher correlation
> rm(list=ls())
> set.seed(123456)
> rho <- 0.75
> 
> accept<-c()
> 
> reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
> reviews <- reviews.all[,1:4]
> 
> # impose desk rejection
> not.desk <- which(reviews[,4] >= 0.5)
> 
> # strong editor system
> decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
> decisions <- decisions.t >= (1-0.09)
> sum(decisions)/length(decisions)
[1] 0.1022945
> 
> review.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)
> 
> lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lwd=2)
> 
> mean(review.avg[which(decisions==T)])
[1] 0.9236401
> 
> 
> 
> legend("topleft", bty="n", cex=0.6, lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black"), legend=c("two subdisciplines, rho = 0.9 or 0.1          ", "all readers correlated at rho = 0.5          ", "all readers correlated at rho = 0.75          ") )
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
 226.25    4.73  231.28 
